In [1]: 
%pylab inline
import scipy
import scipy.fftpack
import scipy.signal

import scipy.io

load_path = "/home/jbm/capture_aprs.raw"
Fc = 144.375e6
Fs = 250e3


pylab.rcParams['figure.figsize'] = (15.0, 6.0)

Populating the interactive namespace from numpy and matplotlib
In [2]: 
samps = [ 512000, 128000, 1024000, 16000 ]

#x = sdr.read_samples(2*sdr.sample_rate)
x = np.fromfile(load_path, dtype="complex64")

figure(figsize=(15,4))
psd(x, Fs=Fs)
title("Raw input")
None
In [3]: 
samps = [ 16000, 128000, 512000, 128000 ]
N = len(samps)

fs, b = scipy.signal.welch(x, fs=Fs, nperseg=1024)
fs = scipy.fftpack.fftshift(fs)
b = scipy.fftpack.fftshift(b)

bl = np.log(b)
dbl = diff(bl)

None
In [4]: 
y = sorted(dbl)

[ y[i*len(y)/100] for i in [10, 25, 50, 75, 90] ]
Out[4]:
[-0.015214484966499242,
 -0.0067623277127637493,
 0.00031054897998572528,
 0.008446068080306901,
 0.016065986410019661]
In [5]: 
il = 5*len(y)/100
ih = 95*len(y)/100
yp = y[il:ih]

my_floor = np.mean(yp) + np.std(yp)*3.0

hits = (np.array(dbl) > my_floor).nonzero()[0]+1
hits = hits[ (hits > len(bl)/10) & (hits < len(bl)*9/10) ]

[ f(yp) for f in [ np.mean, np.std] ]
Out[5]:
[0.00054897979808226087, 0.0095945428733462488]
In [6]: 
# Do some janky RLEish work to find contiguous channels.
# There are better numpy-ish ways to do this, but this is legible.
runs = []

last_start = None

gap_allowed = 2 # Allow some fudge factor here

for (h0,h1) in zip(hits[:-1], hits[1:]):
    if h1 - h0 < gap_allowed:
        if None == last_start:
            last_start = h0
    else:
        if None != last_start:
            runs.append( (last_start, h0) )
            last_start = None
            
            
plot(fs, bl)
for h0,h1 in runs:
    if h1-h0 > 1:
        axvspan(fs[h0-1], fs[h1+1], color="red", alpha=0.2)
xlim(fs[len(bl)/10], fs[9*len(bl)/10])
            
[ ((fs[h0]+fs[h1])/2, (fs[h1]-fs[h0]), fs[h0], fs[h1]) for (h0,h1) in runs if h1-h0 > 1 ]
Out[6]:
[(-4638.671875, 488.28125, -4882.8125, -4394.53125),
 (-2197.265625, 3417.96875, -3906.25, -488.28125),
 (10253.90625, 488.28125, 10009.765625, 10498.046875),
 (12573.2421875, 732.421875, 12207.03125, 12939.453125),
 (14892.578125, 488.28125, 14648.4375, 15136.71875)]
In [7]: 
specgram(x[:Fs/2], Fs=Fs, NFFT=1024, noverlap=512)

for (h0,h1) in runs:
    axhspan(-1+2.0*float(h0)/len(bl), -1+2.0*float(h1)/len(bl), color="grey", alpha=0.5)
None
In [8]: 
def extract_channel(f_sample, f_center, f_bw, x):
    """Extract a channel of width f_bw at f_center, from f_sample
    
    Returns Fs,y, the new sample rate and an array of complex samples
    """
    
    # Cribbing from the GNURadio xlate, we use their block diagram
    
    my_state = {} # cheating backchannel for debug
    
    # Create a LPF for our target bandwidth
    n_taps = 100 # XXX Total random shot in the dark    
    r_cutoff = float(f_bw)/f_sample    
    base_lpf = scipy.signal.remez(n_taps, [0, r_cutoff, r_cutoff*2, 0.5], [1,0])
    
    # Shift this LPF up to our center frequency
    fwT0 = 2.0*np.pi*f_center/f_sample
    
    lpf = base_lpf * np.exp(1.0j*fwT0 * np.arange(0, n_taps))
    
    dec_rate = int(f_sample / (2*f_bw))
    
    x_filtered = scipy.signal.lfilter(lpf, 1.0, x)
    
    dec_is = np.arange(0, len(x_filtered), dec_rate)
    y = x_filtered[dec_is]
    y *= np.exp(-1.0j*fwT0*dec_is)
    
    my_state["n_taps"] = n_taps
    my_state["r_cutoff"] = r_cutoff
    my_state["lpf"] = lpf
    my_state["base_lpf"] = base_lpf
    my_state["fwT0"] = fwT0
    my_state["dec_rate"] = dec_rate
    # my_state["dec_is"] = dec_is
    # my_state["x_filtered"] = x_filtered
    
    return f_sample/dec_rate, y, my_state
In [9]: 
fs_new, y, extract_state = extract_channel(Fs, 15000, 5000, x)

specgram(y[0:100000], Fs=fs_new)

fs_new
Out[9]:
10000.0
In [10]: 
del(extract_state)
In [11]: 
def decode_quad(x, gain):
    xp = x[1:] * np.conj(x[:-1])
    retval = gain * np.arctan2(np.imag(xp), np.real(xp))
    return retval
In [15]: 
gain = Fs / (2.0 * np.pi * fs_new)  # Cribbed from the gnuradio code
w = decode_quad(y, gain)
specgram(w[0:100000], Fs=fs_new) # , NFFT=1024, noverlap=1023)
#axhline(1200)
#axhline(2200)
#xlim(1.75, 2)
fs_new, min(w), max(w)
Out[15]:
(10000.0, -12.499961187501176, 12.499991045392912)
In [13]: 
w_out = w[50:-50] * 500.0
w_out -= np.mean(w_out)
w_out.astype("int16").tofile("foo.raw")
In [14]: 
specgram(w[50:1000], Fs=fs_new); None